import base64
import statistics

import joblib
import numpy as np
import pandas as pd
from scipy import signal
from scipy.signal import welch

ecg_data_all=[]#全局变量，存储所有心电数据
class ModelHandler:
    def __init__(self):
        self.ecg_data = ecg_data_all
        # 绘制心电图
        # print(len(self.ecg_data))
        # plt.figure(figsize=(12, 6))
        # plt.plot(self.ecg_data)
        # plt.title('ECG Signal')
        # plt.xlabel('Sample Number')
        # plt.ylabel('Amplitude')
        # plt.show()
        # print(self.ecg_data )
    def handler_data(self, data):
        if 'ecgData' in data:
            self.handle_ecg_date(base64.b64decode(data['ecgData']), 'CHAN_8', False)
            # print('原始数据：',data['ecgData'])
            # print('\n')
            # print('编码数据：',base64.b64decode(data['ecgData']))

            # print(len(data['ecgData']))
            # print('\n')
        # 模型调用代码补充
        # ……
        # print(len(self.ecg_data))
        # if len(self.ecg_data) > 0:
        #     ecg_data_all.append(self.ecg_data)
        # print(ecg_data_all)
        HRV_data=self.process_mv_data()
        # print("HRV_data",HRV_data)
        if HRV_data:
            # HRV_data =self.process_mv_data()
            print(HRV_data)
            selected_hrv_params = ['HRV_MeanNN','HRV_SDNN','HRV_RMSSD','HRV_NN50','HRV_pNN50','HRV_CVNN','HRV_SDSD','HRV_LFn','HRV_HFn','HRV_LFn/HFn','HRV_SD1','HRV_SD2','HRV_SD1SD2','HRV_S']
            selected_values = {k: v for k, v in HRV_data.items() if k in selected_hrv_params}
            print(selected_values)
            # if selected_values['HRV_LFn/HFn']:
            #     selected_values['HRV_LFn/HFn']=0

        # selected_values = [HRV_data[param][0] for param in selected_hrv_params if param in HRV_data]
            # 加载模型
            df = pd.DataFrame(selected_values)
            # 将无穷大值替换为 0
            df.replace([np.inf, -np.inf], 0, inplace=True)

            # 将 NaN 值替换为 0
            df.fillna(0, inplace=True)
            print(df)
            af_model = joblib.load(open('LsTM-3class-30.pkl','rb'))
            result = af_model.predict_proba(df)
            result=result.flatten()
            print(result)

            return result
        else:
            result=None
            # print(result)
            return result

    # 解析获取到的心电数据
    def handle_ecg_date(self, ecgData, channelType, original):
        dataLen = len(ecgData)
        # print(dataLen)
        index = 0
        if "CHAN_1" == channelType:
            # 单通道模式
            while index < dataLen - 1 * 3 + 2:
                value = 0
                b = ecgData[index]
                value |= (b & 0x000000ff)
                value <<= 8
                index += 1

                b = ecgData[index]
                value |= (b & 0x000000ff)
                value <<= 8
                index += 1

                b = ecgData[index]
                value |= (b & 0x000000ff)
                index += 1

                if original:
                    mv = (value - 0x800000) * 0.476837 / 12 / 1000
                else:
                    mv = (value - 6000) / 820.0
                # print(mv)
                self.ecg_data.append(mv)
        elif "CHAN_3" == channelType:
            # 3通道模式
            while index < dataLen - 3 * 3 + 2:
                ad_value = []
                for i in range(3):
                    value = 0
                    b = ecgData[index]
                    value |= (b & 0x000000ff)
                    value <<= 8
                    index += 1

                    b = ecgData[index]
                    value |= (b & 0x000000ff)
                    value <<= 8
                    index += 1

                    b = ecgData[index]
                    value |= (b & 0x000000ff)
                    index += 1

                    ad_value.append(value)
                if original:
                    mv = (value - 0x800000) * 0.476837 / 12 / 1000
                else:
                    mv = (value - 6000) / 820.0
                # print(mv)
                self.ecg_data.append(mv)
        elif "CHAN_8" == channelType:
            # 8通道模式
            while index < dataLen - 8 * 3 + 2:
                ad_value = []
                for i in range(8):
                    value = 0
                    b = ecgData[index]
                    value |= (b & 0x000000ff)
                    value <<= 8
                    index += 1

                    b = ecgData[index]
                    value |= (b & 0x000000ff)
                    value <<= 8
                    index += 1

                    b = ecgData[index]
                    value |= (b & 0x000000ff)
                    index += 1

                    ad_value.append(value)
                if original:
                    mv = (value - 0x800000) * 0.476837 / 12 / 1000
                else:
                    mv = (value - 6000) / 820.0
                # print(mv)
                # print('输出数据', mv)

                # if mv:
                ecg_data_all.append(mv)

                # print(self.ecg_data)


        # 处理mv数据并计算HRV

    def detect_r_peaks(ecg_signal, fs=500):
        r_peaks = []  # 存储R波位置的索引
        # 这里简化处理，假设直接通过阈值检测R波位置
        threshold = 0.6 * max(ecg_signal)
        for i in range(1, len(ecg_signal) - 1):
            if ecg_signal[i] > threshold and ecg_signal[i] > ecg_signal[i - 1] and ecg_signal[i] > ecg_signal[i + 1]:
                r_peaks.append(i)
        return r_peaks

    def calculate_rr_intervals(r_peaks, fs=500):
        # 计算R-R间期
        rr_intervals = []
        for i in range(1, len(r_peaks)):
            rr_interval = (r_peaks[i] - r_peaks[i - 1]) / fs
            rr_intervals.append(rr_interval)
        return rr_intervals

    def calculate_rr_intervals_from_mv_chunk(mv_chunk, fs=500):
        # mv_chunk是包含心电信号数据点的列表或数组
        r_peaks = ModelHandler.detect_r_peaks(mv_chunk, fs)
        rr_intervals = ModelHandler.calculate_rr_intervals(r_peaks, fs)
        return rr_intervals

    def process_mv_data(self):
        hrv_result = None
        # print(self.ecg_data)
        # self.ecg_data = [data for data in self.ecg_data if data]
        while len(self.ecg_data) >= 0.5 * 60 * 500:
            mv_chunk = self.ecg_data[:int(0.5 * 60 * 500)]
            del self.ecg_data[:int(0.5 * 60 * 500)]  # 删除已处理的数据
            del ecg_data_all[:int(0.5 * 60 * 500)]  # 删除已处理的数据
            # print(mv_chunk)
            # print(len(mv_chunk))
            # # 绘制心电图
            # plt.figure(figsize=(12, 6))
            # plt.plot(mv_chunk)
            # plt.title('ECG Signal')
            # plt.xlabel('Sample Number')
            # plt.ylabel('Amplitude')
            # plt.show()
            # print(mv_chunk.shape)
            # print(type(mv_chunk))
            rr_intervals = ModelHandler.calculate_rr_intervals_from_mv_chunk(mv_chunk)
            print(rr_intervals)
            # 这里调用HRVUtil中的计算函数计算HRV，示例中假设存在HRVUtil类来处理HRV计算
            hrv_util = HRVUtil(rr_intervals)  # 假设存在HRVUtil类来处理HRV计算

            hrv_result = hrv_util.calculate_all_hrv_indices(rr_intervals)
            # 在这里可以将hrv_result存储到数据库或者其他持久化存储中
            # print(hrv_result)
        return hrv_result


class HRVUtil:
    def __init__(self,rr_intervals):
        self.mv_data = rr_intervals  # 存储mv数据的列表

    def calculate_hrv_meanNN(self):
        # 计算HRV_MeanNN，即RR间期的平均值
        meanNN = sum(self.mv_data) / len(self.mv_data)
        return meanNN

    def calculate_hrv_sdnn(self):
        # 计算HRV_SDNN，即RR间期的标准差
        import statistics
        sdnn = statistics.stdev(self.mv_data)
        return sdnn

    # 其他HRV参数的计算函数依次类推

    def calculate_hrv_rmssd(self):
        nn_intervals = np.diff( self.mv_data)
        rmssd = np.sqrt(np.mean(np.square(nn_intervals)))
        return [rmssd]

    def calculate_hrv_sdsd(self):
        sdsd = np.std(np.diff( self.mv_data))
        return [sdsd]

    def calculate_hrv_cvnn(self):
        nn_intervals = np.diff( self.mv_data)
        mean_nn = np.mean(nn_intervals)
        sdnn = np.std(nn_intervals)
        cvnn = sdnn / mean_nn * 100
        return [cvnn]

    def calculate_hrv_cvsd(self):
        sdnn = np.std( self.mv_data)
        mean_rr = np.mean( self.mv_data)
        cvsd = sdnn / mean_rr * 100
        return [cvsd]

    def calculate_hrv_nn50(self):
        # 假设mv_chunk是一个包含RR间期的列表或数组
        rr_intervals = self.mv_data
        # 计算相邻RR间期的差值
        rr_diffs = [rr_intervals[i + 1] - rr_intervals[i] for i in range(len(rr_intervals) - 1)]

        # 统计大于50ms的次数
        nn50_count = sum(1 for diff in rr_diffs if abs(diff) > 50)

        return nn50_count

    def calculate_hrv_pnn50(self):
        nn_intervals = np.diff( self.mv_data)
        pnn50 = np.sum(np.abs(nn_intervals) > 50) / len(nn_intervals) * 100
        return [pnn50]

    def calculate_hrv_pnn20(self):
        nn_intervals = np.diff(self.mv_data)
        pnn20 = np.sum(np.abs(nn_intervals) > 20) / len(nn_intervals) * 100
        return [pnn20]

    def calculate_hrv_lfn(self):

        # 步骤1: 进行傅立叶变换 (FFT)
        fft_result = np.fft.fft(self.mv_data)

        # 步骤2: 计算功率谱密度 (PSD)
        psd = np.abs(fft_result) ** 2

        # 步骤3: 提取LF带宽内的能量
        fs = 500  # 采样率为500 Hz
        f, Pxx = signal.welch(self.mv_data, fs=fs, nperseg=len(self.mv_data))

        # 步骤4: 计算LF功率
        lf_band = (f >= 0.04) & (f <= 0.15)
        lf_power = np.trapz(Pxx[lf_band], f[lf_band])

        return lf_power

    def calculate_hrv_hf(self):
        # Example implementation for HF calculation (frequency domain method)
        f, Pxx = welch(self.mv_data, fs=500, nperseg=len(self.mv_data))
        hf_band = (0.15, 0.4)
        hf_power = np.trapz(Pxx[(f >= hf_band[0]) & (f <= hf_band[1])], f[(f >= hf_band[0]) & (f <= hf_band[1])])
        return [hf_power]

    def calculate_hrv_vhf(self):
        # Example implementation for VHF calculation (frequency domain method)
        f, Pxx = welch(self.mv_data, fs=500, nperseg=len(self.mv_data))
        vhf_band = (0.4, 0.8)
        vhf_power = np.trapz(Pxx[(f >= vhf_band[0]) & (f <= vhf_band[1])], f[(f >= vhf_band[0]) & (f <= vhf_band[1])])
        return [vhf_power]

    def calculate_hrv_hfn(self):
        # Example implementation for HF normalized (HFn) calculation
        f, Pxx = welch(self.mv_data, fs=500, nperseg=len(self.mv_data))
        total_power = np.trapz(Pxx, f)
        hf_band = (0.15, 0.4)
        hf_power = np.trapz(Pxx[(f >= hf_band[0]) & (f <= hf_band[1])], f[(f >= hf_band[0]) & (f <= hf_band[1])])
        hfn = hf_power / total_power * 100
        return [hfn]

    def calculate_hrv_lnhf(self):
        # Example implementation for LnHF calculation
        f, Pxx = welch(self.mv_data, fs=500, nperseg=len(self.mv_data))
        hf_band = (0.15, 0.4)
        hf_power = np.trapz(Pxx[(f >= hf_band[0]) & (f <= hf_band[1])], f[(f >= hf_band[0]) & (f <= hf_band[1])])
        ln_hf = np.log(hf_power)
        return [ln_hf]

    def calculate_hrv_lfn_hfn(self):
        x=self.calculate_hrv_lfn()
        y=self.calculate_hrv_hfn()
        return x/y

    def calculate_hrv_sd1(self):
        # SD1
        sd1 = np.std(np.diff(self.mv_data)) / np.sqrt(2)
        return [sd1]

    def calculate_hrv_sd2(self):
        # SD2
        mv_data_array = np.array(self.mv_data)
        sd2 = np.std(mv_data_array[:-1] - mv_data_array[1:]) / np.sqrt(2)
        return [sd2]

    def calculate_hrv_sd1sd2(self):
        # 计算SD1/SD2比值
        mv_data_array = np.array(self.mv_data)
        sd1 = np.std(np.diff(self.mv_data)) / np.sqrt(2)
        sd2 = np.std(mv_data_array[:-1] - mv_data_array[1:]) / np.sqrt(2)
        sd1_sd2_ratio = sd1 / sd2
        return [sd1_sd2_ratio]

    def calculate_hrv_s(self):
        # 香农熵（Shannon Entropy，S）
        p = np.histogram(self.mv_data, bins='auto')[0] / len(self.mv_data)
        s = -np.sum(p * np.log2(p + 1e-12))
        return [s]

    def calculate_hrv_apen(self):
        # 近似熵（Approximate Entropy，ApEn）
        m = 2  # Embedding dimension
        r = 0.2 * np.std(self.mv_data)  # Tolerance threshold
        N = len(self.mv_data)

        def _phi(m):
            x = np.array([self.mv_data[i:i + m] for i in range(N - m + 1)])
            return np.sum(np.max(np.abs(x[:, np.newaxis] - x[np.newaxis, :]), axis=2) <= r) / (N - m + 1)

        apen = -np.log(_phi(m) / _phi(m + 1))
        return [apen]

    def calculate_meanHR(self):
        # 计算心率的均值
        meanHR = 60 / self.calculate_hrv_meanNN()
        return meanHR

    def calculate_sdHR(self):
        # 计算心率的标准差
        rr_intervals = [60 / rr for rr in self.mv_data]
        sdHR = statistics.stdev(rr_intervals)
        return sdHR

    def calculate_sampen(self, m, r):
        # 计算样本熵（Sample Entropy, Sampen）
        def _phi(m):
            x = np.array([self.mv_data[i:i + m] for i in range(len(self.mv_data) - m + 1)])
            C = np.sum([np.sum(np.abs(x - x[j]).max(axis=1) <= r) - 1 for j in range(len(x))]) / (len(x) - 1)
            return C

        return -np.log(_phi(m + 1) / _phi(m))

    def calculate_alpha1(self):
        # 计算去趋势波动分析的短程标度指数（alpha1）
        def dfa(data):
            N = len(data)
            all_n = range(4, int(np.floor(N / 4.0)))
            F = []
            for n in all_n:
                rms = []
                for i in range(0, N - n, n):
                    c = np.polyfit(range(n), data[i:i + n], 1)
                    fit = np.polyval(c, range(n))
                    rms.append(np.sqrt(np.mean((data[i:i + n] - fit) ** 2)))
                F.append(np.mean(rms))
            alphas = np.polyfit(np.log(all_n), np.log(F), 1)
            return alphas[0]

        return dfa(self.mv_data)

    def calculate_all_hrv_indices(self,rr_intervals):
        # 计算所有的HRV参数
        hrv_indices = {
            'HRV_MeanNN': self.calculate_hrv_meanNN(),
            'HRV_SDNN': self.calculate_hrv_sdnn(),
            'HRV_RMSSD': self.calculate_hrv_rmssd(),
            'HRV_SDSD': self.calculate_hrv_sdsd(),
            'HRV_CVNN': self.calculate_hrv_cvnn(),
            'HRV_CVSD': self.calculate_hrv_cvsd(),
            'HRV_NN50': self.calculate_hrv_nn50(),
            'HRV_pNN50': self.calculate_hrv_pnn50(),
            'HRV_pNN20': self.calculate_hrv_pnn20(),
            'HRV_HF': self.calculate_hrv_hf(),
            'HRV_LFn': self.calculate_hrv_lfn(),
            'HRV_VHF': self.calculate_hrv_vhf(),
            'HRV_HFn': self.calculate_hrv_hfn(),
            'HRV_LnHF': self.calculate_hrv_lnhf(),
            'HRV_LFn/HFn': self.calculate_hrv_lfn_hfn(),
            'HRV_SD1': self.calculate_hrv_sd1(),
            'HRV_SD2': self.calculate_hrv_sd2(),
            'HRV_SD1SD2': self.calculate_hrv_sd1sd2(),
            'HRV_S': self.calculate_hrv_s(),
            'HRV_ApEn': self.calculate_hrv_apen(),
            'HRV_meanHR': self.calculate_meanHR(),
            'HRV_sdHR': self.calculate_sdHR(),
            'HRV_Sampen': self.calculate_sampen(m=2, r=0.2 * np.std(rr_intervals)),
            'HRV_alpha1': self.calculate_alpha1()
            # 其他HRV参数依次添加
        }
        return hrv_indices